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We study critical hysteresis in the random-field Ising model (RFIM) on a two-dimensional periodic 
lattice with a variable coordination number Zeff in the range 3 < Zeff < 6. We find that the model 
supports critical behavior in the range 4 < z^ff < 6, but the critical exponents are independent of 
Zef f ■ The result is discussed in the context of the universality of nonequilibrium critical phenomena 
and extant results in the field. 


I. INTRODUCTION 

Disordered systems have been studied extensively over the last few decades. A point of special interest is the effect of 
quenched disorder in the system. The quenched disorder endows the system with a large number of metastable states 
that are surrounded by high energy barriers. This alters how the system relaxes and how it responds to an external 
driving force. The response becomes sporadic and jerky even if the driving force increases slowly and smoothly. The 
reason is that the system is unable to move from one metastable state to another unless the applied force has reached 
the necessary level for crossing the barrier between the two metastable states. The passage from one state to another 
involves an avalanche (a quick succession of restructuring events) in the internal structure of the system. A large 
variety of systems exhibit this avalanche dynamics. Examples include martensitic transformations in many metals and 
alloys as they are cooled from a bcc structure to a closed packed structure [l|, earthquakes caused by the movement 
of tectonic plates 0, re-arrangement of domains in a ferromagnet, and motion in granular materials Q. The size 
of an avalanche may vary over a wide range; it may be microscopic, large, or critical in the sense of a diverging size. 
We are particularly interested in this paper in nonequilibrium critical phenomena which are caused by a diverging 
avalanche. This is akin to critical phenomena in systems in equilibrium that are caused by a diverging correlation 
length. There is a good deal of experimental support largely from Barkhausen noise experiments with a wide variety 
of materials 0,0 that nonequilibrium critical phenomena have much in common with their equilibrium counterpart 
such as power laws spanning many decades and universal critical exponents. 

Apparently there is a significant difference between universality of equilibrium critical phenomena and its coun¬ 
terpart in the avalanche-driven behavior. In equilibrium, the dimensionality d of the system is a key player that 
determines the universality class of critical behavior. Short-range structure of the system is irrelevant for this pur¬ 
pose. The rationale is that the critical behavior is caused by spontaneous fluctuations in the system whose size diverges 
at the critical point. It should not be influenced by short-range details of the system. Thus the critical phenomena 
on a sc, bcc, or fee structure should be the same. This idea is well tested theoretically as well as experimentally in 
the equilibrium case and it may be expected to hold in the nonequilibrium case as well. However we shall see in the 
following that it is not the case. Extensive study of nonequilibrium critical behavior has been carried out in the frame¬ 
work of hysteresis in the random-field Ising model (RFIM) on a sc lattice with on-site quenched random-field having a 
Gaussian distribution with mean value zero and standard deviation a 00. The basic quantity studied is m{h, a), the 
magnetization per site in the system at zero temperature, as the applied field h is ramped up slowly from h = —oo to 
h = oo. Extensive numerical studies on the sc lattice [7 h 10I| reveal a critical value of cr = ctc ~ 2.16 such that m(h, a) 
is macroscopically continuous if cr > CTc but has a jump discontinuity if cr < CTc- The size of the jump discontinuity 
reduces as cr CTc and the jump occurs at larger values of the applied field h. The point (ctc = 2.16, he = 1.435) is 
a nonequilibrium critical point marked by scale invariant phenomena similar to the equilibrium critical phenomena. 
The model has been studied on a square lattice as well pTl - ll3l| . Initial results on the square lattice were inconclusive 
but the most recent studies [10 show the existence of a critical point. These results suggest that 2d may be the lower 
critical dimension for the nonequilibrium critical behavior. 

The problem of hysteresis in the nonequilibrium random-field Isii^ model at zero temperature can be solved 
analytically on a Bethe lattice of an arbitrary coordination number z [l4| . The exact solution brings out a surprising 
fact that the nonequilibrium critical point (tTc, he) does not exist if 2 ; < 4. Numerical studies of the model on periodic 
lattices suggest a similar result i.e. the absence of avalanche-driven criticality if 2 ; < 4, irrespective of the spatial 
dimensionality d of the lattice [T0. For example, consider two different 2d periodic lattices, (i) the honeycomb lattice 
with z = 3, and (ii) the triangular lattice with z = 6. A nonequilibrium critical point is absent on the honeycomb 
lattice [l0 but present on the triangular lattice [10. The pattern of extant results suggests that a lower critical 
coordination number (ze = 4) has greater significance than the idea of a lower critical dimension for nonequilibrium 
critical phenomena. In order to examine this point further, we study the problem on a randomly diluted triangular 
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lattice. The triangular lattice comprises three equivalent inter-penetrating sublattices A, B, and C. We randomly 
decimate the sites on one of these sublattices, say the sublattice C. If c is the concentration, i.e. the fraction of sites 
present on the sublattice C, then the limit c = 0 corresponds to the two sub-lattices A and B making up a honeycomb 
lattice with coordination number equal to 3. For other values of the dilution parameter c, sites on the lattices A or B 
have an average coordination number Zeff{c) = 3(1 + c). We look for the presence of a critical point in the hysterestic 
response of the system for various values of c and find that the critical point disappears if c < 0.33, or equivalently if 
the effective coordination number Zeff{c) is less than 4. We also examine the dependence on c of the critical exponent 
v that characterizes the divergence of the largest avalanche on the lattice. We conclude that v is independent of c in 
the range (0.33 < c < 1) within numerical errors. 

The paper is organized as follows. In section-II we describe the RFIM and its dynamics very briefly for the sake 
of completeness and to set up our notation. Section III presents numerical results for the magnetization curve on the 
lower branch of the hysteresis loop for a selected value of dilution (c = 0.90), and selected values of a characterizing 
the random-field distribution {a = 0.95, 1.01, and 1.05). This section emphasizes the key idea for determination of 
the critical point. Section IV contains the body of our numerical results for a systematically diluted lattice and finite 
size analysis. Section V contains a brief discussion of the main results of our study. 


II. THE MODEL AND QUALITATIVE BEHAVIOR 

The RFIM with nearest neighbor ferromagnetic interaction J is characterized by the Hamiltonian, 


H— j'Y2,CiCjSiSj '^hiCiSi h'Y2,CiSi (1) 

i i 

Here Ci is a quenched binary variable; = 1 if site i oi a. L x L triangular lattice shown in figure (1) is occupied by 
an Ising spin Si = ±1, and Ci = 0 otherwise. Each site has a quenched random-field hi with a Gaussian distribution of 
mean value zero and standard deviation a. The system is placed in a uniform external field h which is varied slowly 
from h = —oo to h = oo. 

The discrete time, single spin flip serial dynamics of the zero-temperature nonequilibrium RFIM is specified by the 
equation El, 


Si(t -|- 1) = sign h = J cjSj + hi + h (2) 

3 

Figure (1) shows three inter-penetrating sub-lattices A, B and C that comprise a triangular lattice. We randomly 
dilute the sub-lattice C. Thus only a fraction c of the sites on the sublattice C are occupied by Ising spins. It is easily 
seen that in the limit c = 0, the triangular lattice reduces to the honeycomb lattice and the coordination number 
of the lattice reduces from 6 to 3. For other values of c in the range 0 < c < I, the lattice is characterized by an 
inhomogeneous coordination number. An occupied site on the sublattice C always has 6 nearest neighbors; 3 on 
sublattice A and 3 on sublattice B. A site on sublattice A (B) has 3-|-3c nearest neighbors; 3 on sublattice B (A) 
and 3c (on average) on sublattice C. Thus the coordination number for an occupied site on C is 6, and the effective 
coordination number on A or B is given by Ze//(c) = 3(l-|-c). The average coordination number when A, B, and C 
are all taken together is Zav = 6(1 -I- 2c)/(2-|-c). In the following we study the magnetization and the character of the 
avalanches on the diluted triangular lattice for different values of dilution c. The magnetization per site m(/i, a, c; t) 
on the lattice is given by 


m{h, a, c; t) 


1 

Lx L 


LxL 

CiSi 


I 


(3) 


We keep h fixed and iterate the single spin flip dynamics till it reaches a fixed point i.e. a t-independent magne¬ 
tization. We start with m{h = —oo, cr, c) = —(2 -|- c)/3 and raise h till some spin becomes unstable and flips up. It 
may cause some of its neighbors to flip up as well. The dynamics is iterated till a stable configuration is reached. 
The magnetization of the stable configuration is calculated and the process is repeated by raising h to the next higher 
value that makes a new site unstable. The number of spins that flip up in going from one fixed point to the next is 
the size of the avalanche. 
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The size of an avalanche depends on three quantities; the value of the dilution parameter c, the standard deviation a 
of the quenched random field, and the applied field h. For c = 1 there is a critical value of ct = crd^) that separates 
two different behaviors that are easy to understand qualitatively. For cr >> crc(l), the distribution of random-field 
is very wide. Therefore the spins tend to flip up independently of each other and avalanches are relatively small at 
any applied field h. On the other hand, for cr << crc(l) the distribution of on-site random fields is very narrow. In 
this case, a single spin flip that changes the local field at each nearest neighbor by an amount 2\J\ causes a spanning 
avalanche across the system if the applied field is greater than some threshold. This results in a large change in the 
magnetization of the system, i.e. a first-order jump in the magnetization at some critical value of the applied field. 
The size of the jump decreases with increasing a and reduces to zero at the critical point a = crc(l) and h = hc{l). 
We call it a nonequilibrium critical point because the hysteretic susceptibility of the system diverges at this point. 
We find that a similar critical point is present on the diluted lattice in the restricted range 0.33 < c < 1. The value 
of crc(c) reduces with decreasing c and drops discontinuously to zero at c = 0.33 approximately. As mentioned earlier 
c = 1/3 corresponds to Zg// = 4. Thus a critical point occurs on the diluted lattice only if Zeff > 4. This is interesting 
in the context of a similar conclusion reached in the study of the Bethe lattice as well as some other periodic lattices. 

The central object of the present study is the determination of (Jc(c). It seems nearly impossible to determine it 
analytically. This is not surprising given that exact solutions of Ising models, particularly with quenched disorder, 
are rare. We have to resort to numerical simulations of the model. Simulations too have difficulties of choosing a 
good algorithm that suits available computing resources. For each c, we first determine a range (Tmm < <J < Umax 
that may contain Ucic). This is done by locating Umax where no avalanche on the hysteresis loop is macroscopic, and 
a Umin where there is clearly a macroscopic jump in the magnetization. The next step is to study the system for 
different values of a in the range [umimUmax] at suitably chosen intervals 5<j with a view to narrow this range and 
locate the critical point that separates the two behaviors. There are two difficulties. For each value of cr, the avalanche 
distribution is determined by taking an average over a large number of configurations. It is a cpu intensive exercise and 
a compromise is necessary in choosing an optimal 5a. The size of the step 5a is one factor that contributes to the error 
bars in the results. The other difficulty is that fluctuations increase as we approach the vicinity of the critical point. 
The avalanches that were microscopic for a » a^ grow in size and compete with the avalanches associated with a 
jump in magnetization even as the jump approaches zero. We need a method to distinguish the infinite avalanche 
associated with a first order jump in magnetization from an infinite avalanche associated with a critical point. These 
two categories of avalanches are distinguished from each other by looking at their probability distributions. This 
is again a cpu intensive process. Qualitatively the largest avalanche associated with a magnetization jump scales 
linearly with the system size L while the largest avalanche associated with a critical point scales as . The shapes 
of the avalanche distribution functions for a < a^ and cr « CTc are different from each other This difference is 

exploited to pin down Uc within error bars as described in the following. 


III. SIMULATIONS 

The shape of magnetization curve m{h, a, c) in a changing field h depends on the starting state of the system in 
addition to h, cr, and c. We start from m{h = —oo, cr, c) = —(2 -|- c)/3, and raise the field slowly to h = oo. Figure (2) 
shows m{h, a, c) on a 6000 x 6000 lattice for c = 0.90, and three different values of cr = 0.95 (blue triangles), a = 1.01 
(black circles), and a = 1.05 (red squares) in the range 1.5 < h < 1.7. It illustrates three categories of behavior: 
(i) at smaller values of a (blue triangles) m{h, a, c) has a jump discontinuity, (ii) at larger values of a (red squares) 
m{h, cr, c) looks apparently smooth if we allow for fluctuations due to finite size effects that are naturally present in 
any numerical simulation, (iii) there is an intermediate region (black circles) where it is relatively difficult to decide 
if the curve is smooth or has a discontinuity. Somewhere in this region lies a critical value a^ with the corresponding 
magnetization curve being theoretically smooth but containing a point of inflexion at h = he. The inflexion point is 
the critical point where the fluctuations in the system become anomalously large and consequently the susceptibility 
of the system diverges. The difficulty is that it is not easy to identify the inflexion point in simulations. Fluctuations 
also increase with increasing size of the system, and simulated trajectories remain qualitatively similar to the ones 
shown in figure (2). Thus it is unreasonable to expect that simulation with a larger system (which would anyway 
require an unreasonably long computer time) would make it any easier to locate the inflexion point by merely looking 
at the curve. Rather it has to be inferred indirectly from the analysis of fluctuations in its vicinity. The fluctuations 
are anomalously large not only at the critical point {hc,ac) but also at the discontinuity in mjh, a, c) at a < Uc and 
h < he. This presents an additional difficulty in determining the critical point. Our approach [l6l. Il8l| is based on the 
idea that the character of fluctuations at the critical point is different from the character of fluctuations at a first order 
discontinuity. The macroscopic discontinuity in magnetization constitutes a large avalanche but avalanches elsewhere 
on the magnetization trajectory are exponentially small. Thus on a logarithmic scale, the probability P{s,a,c) of 
an avalanche of size s has two parts, a part that decreases linearly with s and another a delta function peak at 
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s ~ Smax] Smax IS of the Order of the system size but decreases as a approaches CTc from below. Also the distribution 
of avalanches is asymmetric on the two sides of the discontinuity. Avalanches tend to be larger at the onset of the 
discontinuity. After the discontinuity has occurred, most of the spins in the system have turned up and the number 
of potential sites that can turn up reduces drastically. Thus avalanches immediately after a jump in magnetization 
tend to be smaller in comparison to those just before it. In contrast to this, avalanches on both sides of a critical 
point are similar to each other. The maximum size of a critical avalanche scales as the square root of the system 
size while it scales linearly with the system size at a discontinuity. However in spite of these distinguishing features 
between a first-order and a second order transition, the issue still remains difficult to decide because as cr ^ (Tc from 
below, the size of the first order jump tends to zero and the Delta function peak at s = Smax tends to vanish. In the 
absence of a more efficient method, we use the same (rather laborious) method to locate critical points on the diluted 
triangular lattice as was used on the undiluted lattice (c = 1). In the following we discuss the case c = 0.90 in detail, 
and present the results for other values of dilution in the form of a table and graphs. 

As discussed in reference 0, we count the occurrence of all avalanches of size s as the system is driven from 
h = — oo to oo. Let P{s,a,c) be the probability of an avalanche of size s anywhere on the m{h,a,c) curve in an 
increasing applied field {h = —oo to h = oo) on a triangular lattice whose sites on sublattice C are occupied with 
probability c. As mentioned earlier, logP(s, cr, c) will have a linearly decreasing part in the range I < s < Smax if 
m{h, a, c) has a discontinuity. This is born out by the red (filled circles) curve in figure (3) which shows the raw data 
for a = l.I on a 240 x 240 lattice for c = 0.90 and 50000 independent realizations of the random-field distribution. 
The green (filled squares) curve in the same hgure shows the data for a = 1.295 which we estimate to be the critical 
value. Simulations were performed for several closely spaced values of a. However we show the data for only three 
values of a so as not to crowd figure (3). The data depicted in green is the closest to a linear decrease. The Delta 
function peak has vanished signifying that the jump in the magnetization has approached zero. Data for cr = 1.5 
is shown in blue (filled triangles). We conclude that the value of a corresponding to the blue must be greater than 
(Tc due to two reasons: (i) the largest avalanche Smax with a non-zero probability of occurrence is far less than the 
system’s size, and (ii) P{s,a,c) does not approach to zero linearly at Smax but rather bends down to it. 

Figure (4) shows similar data as figure (3) but in a binned form. Binning reduces the scatter of the data and we 
are able to show the data for a larger set of closely spaced cr values without crowding the figure too much. The raw 
data for each value of cr has been binned in 50 linear bins. Here the values of a are very close to CTc so that the largest 
avalanche in each case is of the same order of magnitude. What distinguishes different values of a is that curves 
that bend up near the largest avalanche indicate a tendency to form a Delta function peak. These indicate that the 
corresponding cr is smaller than Uc. Similarly the curves that bend down near the largest avalanche indicate that the 
corresponding values of cr are larger than The curve closest to a straight line (blue circles) belongs to (Jc = 1.295. 
Before closing this section, two comments on the binning procedure may be in order. We have also used logarithmic 
binning which is normally the preferred binning procedure when the distributions have a fat tail as in our case. This 
does not change the results presented here. The main purpose of figure (4) is to distinguish between curves that bend 
up near the end from those that bend down. This is easier for the eye if linear binning is used due to higher density of 
data points in the region. The second point is that the last bin in each case should be ignored due to lack of sufficient 
number of data points in this bin. 


IV. FINITE SIZE EFFECTS AND RESULTS 


Using the procedure outlined above, we have determined the critical value ac{L,c) on diluted L x L lattices for 
various values of L in the range 100 to 600. These are rather small sizes compared with the ones used in the study of the 
undiluted problem on a square lattice BUI- However, such small sizes were found to be adequate in reference [l^ 
to demonstrate the presence of critical hysteresis on a triangular lattice. In the present study as well, these appear 
to be adequate to examine the effect of dilution parameter c and the effective coordination number Zeff on critical 
hysteresis in the system. Our results are presented in Table I. According to the scaling hypothesis, the correlation 
length in the vicinity of the critical point scales as ~ [o'c(c) — as cr —crc(c) from below. The argument on zc(c) 

indicates that we are open to the possibility that the the critical exponent ly may depend on the amount of dilution 
on the lattice. However, as we shall see in the following, this turns out not to be true within the error bars of our 
analysis. In the present problem, the correlation length is measured by the largest distance that an avalanche travels 
from its point of origin. On a finite lattice, the farthest distance an avalanche can travel is limited by the size of 
the lattice. Thus L ~ [<7c{L,c) — where ac{L,c) is a lattice-dependent critical value of crc(c). Allowing for a 

constant of proportionality, we may write [l^ 




ac{L,c) - acjc) 
(Jc{c) 


or 


~r~\ 

iy{c) 


log 


10 


=(i,c) 


crc(c) 



(4) 
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Here the factor l/crc(c) appearing immediately after the first equality sign is a constant of proportionality. A more 
general constant of proportionality 1/a is considered in equation (5). 


TABLE I: ac{L,c) 


L 

c = 1.00 

c = 0.90 

c = 0.80 

c = 0.70 

c = 0.60 

c = 0.50 

c = 0.40 

c = 0.34 

99 

1.63 ±0.01 

1.47±0.01 

1.305±0.005 

1.125 ±0.005 

0.885±0.005 

0.665±0.005 

0.51±0.01 

0.485±0.005 

120 

1.59 ±0.01 

1.425±0.005 

1.255±0.005 

1.065 ±0.005 

0.83±0.01 

0.62±0.01 

0.485±0.005 

0.46±0.01 

141 

1.56 ±0.01 

1.39±0.01 

1.215±0.005 

1.02 ±0.01 

0.785±0.005 

0.585±0.005 

0.465±0.005 

0.44±0.01 

168 

1.525 ±0.005 

1.355±0.005 

1.175±0.005 

0.98 ±0.01 

0.745±0.005 

0.555±0.005 

0.445±0.005 

0.425±0.005 

198 

1.50±0.01 

1.325±0.005 

1.14±0.01 

0.94±0.01 

0.71±0.01 

0.525±0.005 

0.425±0.005 

0.405±0.005 

240 

1.47±0.01 

1.295 ±0.005 

1.105±0.005 

0.90±0.01 

0.675±0.005 

0.495±0.005 

0.41±0.01 

0.39±0.01 

300 

1.44±0.01 

1.26±0.01 

1.065±0.005 

0.86±0.01 

0.635±0.005 

0.465±0.005 

0.39±0.01 

0.375±0.005 

360 

1.42±0.01 

1.235±0.005 

1.04±0.01 

0.83±0.01 

0.605±0.005 

0.445±0.005 

0.375±0.005 

0.365±0.005 

390 

1.41±0.01 

1.225±0.005 

1.03±0.01 

0.815±0.005 

0.595±0.005 

0.435±0.005 

0.37±0.01 

0.36±0.01 

480 

1.39±0.01 

1.20±0.01 

1.00±0.01 

0.785±0.005 

0.565±0.005 

0.415±0.005 

0.355±0.005 

0.345±0.005 

600 

1.37±0.01 

1.18±0.01 

0.975±0.005 

0.76±0.01 

0.54±0.01 

0.395 ±0.005 

0.345±0.005 

0.335±0.005 


For each value of dilution c, the data in Table I is used to plot — log^g L vs. \ogiQ[(Tc{L,c)/<7c{c) — 1] for different 
values of parameter (Tc(c). The shape of the plot changes from concave up to a straight line, and then to concave 
down as CTc(c) is increased from zero. We search for the best value of the CTc(c) that produces a straight line. However, 
the transition from concave up to concave down is not sharp and spreads over a broad range. This creates a relatively 
large uncertainty in the value of ac(c) that produces the best fit, and also an uncertainty in the slope of the line. For 
example, the best fit to a straight line for the data in the second column (c = 0.90) is obtained for CTc = 1.00 ± 0.04. 
This is shown in figure (5). The slope of the straight line is equal to l/i^(c) = 0.54 ± 0.07, or i/(c) = 1.85 ± 0.26. 
Similarly, we determine z/(c) for other values of dilution c. Our results for (Tc(c) and i/(c) for different values of dilution 
on the triangular lattice are summarized in the following table. 


TABLE II: (Tc(c) and i/(c) for different values of dilution c 



c= 1.00 

c = 0.90 

c = 0.80 

c = 0.70 

c = 0.60 

c = 0.50 

c = 0.40 

c = 0.34 

CTc 

1.22±0.04 

1.00±0.04 

0.77±0.05 

p 

O 

HO 

O 

HO 

p 

O 

CO 

CO 

o 

0.26±0.04 

0.25±0.03 

0.25±0.03 

t/(c) 

1.78±0.29 

1.85±0.26 

1.87±0.29 

1.83±0.26 

1.85±0.28 

1.64±0.28 

1.76±0.35 

1.79±0.33 


Our simulations reveal a critical value Cc of the dilution parameter; Cc = 0.33 approximately. At c = Cc, crc(c) drops 
to zero abruptly, and remains zero for c < Cc as shown in Figure (6). There is no first order jump in m{h,a,c) for 
c < Cc- We have verified this directly from the simulations as well as inferred it from the following. For c < 0.33, 
the plot — log^gL vs. log;^Q[crc(7^, c)/crc(c) — 1] is never concave up for any crc(c) > 0. This indicates the absence of a 
disorder driven jump in the magnetization and therefore the disappearance of a critical point if c < 0.33. To further 
validate this point, we also tried as in ref [l^ . the scaling form 

.-1 _ cre(L,c) -crc(c) 


a 


(5) 
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where a is an arbitrary parameter. The role of a is to shift the curve along the y-axis without changing its shape. 
If we set (Tc(c) = 0 in the above equation, we find that the plot is always concave up for any value of c > c^. This 
means that for c > Cc, the system possesses a critical point with CTc > 0. On the other hand, for c < Cc, the curve is a 
straight line indicating the absence of a critical point. 

The variation of adc) vs. c has been plotted in figures (6), and iy{c) vs c in figure (7). We may draw the following 
conclusions from these figures. As the lattice is diluted increasingly, it continues to support a critical point but 
the value of adc) decreases. This is intuitively reasonable because the random dilution of the lattice amounts to a 
positional disorder in the system that supplements to the disorder due to the random-field. Thus the critical point on 
the diluted lattice corresponds to a narrower distribution of the random-field as compared with the undiluted case. 
We also note that the critical dilution Cc = 0.33 corresponds to z^ff « 4. 

Figure (7) shows that the exponent i/(c) is independent of c within the error bars. We may expect this universality 
to hold for other exponents as well because of the relationship between different exponents. It is unlikely that only 
one exponent in an equation is universal while others are not. We have tried to extract a bit more information from 
our numerical data regarding other exponents. We have tried the collapse of the integrated avalanche size distribution 
with the following scaling form [l3, [HI : 


P{S, a) - (6) 

Here a is the exponent describing the largest cut-off avalanche i.e. Smax where r = , r is the 

avalanche size exponent, /3 gives the scaling of the change in magnetisation due to spanning avalanche. Am ^ |r|^, 
at the critical field. The exponent 8 describes the scaling of the reduced magnetisation with reduced magnetic field, 
m ~ at Oc- Thus the product P{S,a) x is a function of a single variable S'|r|^/“. Fig(8) shows P(S'|r|^/“) 

vs. for c = 0.80; L = 168, 240 and 390; and different values of a ranging from 1.15 to 1.40. The curves collapse 

on each other reasonably well if we choose r -|- af3S = 2.05 and a = 0.12. 

We have also examined how the exponents extracted from the collapsing curves depend on L and |r|. For this 
purpose, we fix c and L (say c = 0.80, and L = 390), and choose three closely spaced values of cr (say a = 1.05,1.06, 
and 1.07). The three values of a correspond to three closely spaced values of |r|. Let \r\avg denotes the average of 
the three |r| values. We search for the values of r -|- afiS and a that produce the best collapse of the three curves 
associated with |r|a„g. Next we choose several different triplets of closely spaced cr values and determine the exponents 
over a range of values of \r\avg- This exercise is repeated for systems of different sizes (L = 168,240,390,600), and 
for a range of values of c in the range c > 0.33. The results for c = 0.80 are shown in figures (9) and (10) along with 
the range of variations (error bars) in the values of the exponents. These figures also show the extrapolated values 
of the exponents in the limit \r\avg —?> 0 or L —>■ oo. Thus we obtain r + aPS = 2.05 ± 0.05 and a = 0.11 ± 0.02 for 
c = 0.80. We have checked that the values of the exponents for other values of dilution c lie in the same range as for 
c = 0.80. This leads us to conclude that the values of the exponents r + aPS and a are independent of c if c > 0.33. 


V. DISCUSSION 

We have investigated critical hysteresis on a two dimensional lattice with a variable coordination number. We had 
to study the problem numerically because exact analytical solutions of this problem are not possible. Simulations of 
the model suggest that the three parameters of the model c, h, and a are characterized by critical values Cc, he and 
Gc respectively. If c > Cc, the critical point [he-, ad is marked by a diverging avalanche and a critical exponent ;z(c). 
Within numerical errors, we find that the exponent v does not depend on c. The critical value Cc is approximately 
equal to Cc = 1/3 which corresponds to Zeff = 4 and Zav = 30/7. We note that each site on A or B sublattice has three 
nearest neighbor sites on the C sublattice which are occupied independently with probability c. Therefore if c < 1/3, 
the probability of a spanning path on A-l-B through occupied sites on C goes to zero. In this case the C sublattice 
does not contribute to a diverging correlation length. The cooperative behavior of the system is qualitatively the 
same as on the honeycomb lattice comprising A and B, and the relevant parameter is Zeff rather than Zav 

We note that c = 1/3 corresponds to an effective coordination number Zeff = 4. Thus our results suggest that 
critical behavior disappears when the coordination number of the lattice drops below four. A similar result holds for 
the Bethe lattice [l^ and also some periodic lattices [l^ . On the Bethe lattice, the problem can be solved analytically 
and therefore the mathematical reason for the absence of critical behavior on Bethe lattice of coordination number 
less than four is understood. However the physical reason for this behavior on periodic lattices is not well understood. 
It is curious that the lower critical coordination number is also equal to four on the diluted triangular lattice. We 
could not have anticipated this result beforehand. Our numerical results only suggest Zeff to be approximately four 
within errors but the probabilistic argument mentioned above indicates that it may be exactly equal to four. Thus 
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all the extant results indicate that nonequilibrium critical behavior occurs only on lattices with coordination number 
equal to four or more irrespective of the dimensionality of the lattice. 

Before closing, we wish to comment on the small difference between results presented here for the case c = 1, and 
those in reference (16). We find CTc = 1.22±0.04 and v = 1.78±0.29 (for c=1.00) while the values reported in reference 
(16) are CTc = 1.27 and v = 1.6 ± 0.2. The results are consistent with each other within error bars. As authors of 
reference (16) are also coauthors of the present paper, we can point out the reason for the difference in the two sets 
of results. We estimate CTc as the best value of this parameter that results in a straight line when the left-hand-side of 
equation (4) is plotted against its right-hand-side. As explained in the previous section, the range of values of crc(c) 
that produces an apparent straight line is rather broad. The results in reference (16) rely on a visual scrutiny of the 
plots and choosing the fit that looks best to the eye. This is of course susceptible to human error. In the present 
study, we have used the linear least squares fitting technique. It has the convenience of a mechanical method to handle 
the data contained in Table I but it also has the disadvantage that outlying points have a disproportionate effect on 
the fit. This contributes to the uncertainty in (Tc(c) and v which is determined by the slope of the straight line. The 
uncertainties in our results are relatively large for the effort put in this study. This seems unavoidable with the present 
method. The remark on the closeness of the exponent v on the triangular lattice [iBl and the simple cubic lattice 
M should also be taken in the same vein. At a qualitative level, it does suggest that the avalanche-driven exponent 
may depend on the coordination number rather than the dimensionality of the lattice. However, it is difficult to reach 
a stronger conclusion with the error bars in our analysis. We hope the present study will motivate further studies of 
the issues raised here with new and improved techniques. 
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FIG. 1: Triangular lattice with sublattices j4(black circles), B (shaded circles) and (^(white circles). C is randomly diluted. 
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FIG. 2: (color online) Magnetization curves on a 6000 x 6000 diluted triangular lattice with c = 0.90 and for cr = 0.95(blue 
triangles), a = 1.01(black circles), and a — 1.05(red squares). 
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FIG. 3: (color online) Distribution of avalanches P{s,a,c) on a. L x L lattice for L — 240 and c = 0.90; a — 1.10 (red circles 
with a peak at s « L x L), a — 1.295 (green squares, note the peak at s ~ L x L has vanished and the plot is almost linear) 
and cr = 1.50 (blue triangles, the bending of the curve indicates the largest avalanche is smaller than L x L). 



FIG. 4: (color online) P{s, a, c) for c = 0.90 and L — 240 using linear binning; a = 1.25(red squares), a — 1.27(green triangles), 
cr = 1.295(blue circles), a = 1.31(pink inverted triangles), and a = 1.35(brown diamonds). Figure shows ac{L,c) « 1.295 
because the corresponding curve is nearly linear; the curves for cr < crc(L, c) tend to peak at the largest avalanche while those 
for cr > ac{L, c) tend to bend down. 
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FIG. 5: (color online) Plot of logio[2i^^(^ — 1] vs. -logic L for c = 0.90. The best fit to a straight line is obtained for 
(Tc(c) = 1.00 ± 0.04. The slope of the straight line yields v[c) = 1.85 ± 0.26. 



FIG. 6: (color online) Variation of crc(c) with c. The figure shows the absence of critical hysteresis if c < 0.33. 
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FIG. 7: (color online) Variation of the exponent z^(c) with c; z^(c) « 1.8 ± 0.3 for c > 0.33 
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FIG. 8: (color online) Scaling collapse of the integrated avalanche size distribution for c = 0.80; L = 168, 240, and 390 for seven 
values of o chosen from the range 1.15 — 1.40. The best collapse is obtained for r + a/35 = 2.05 and a = 0.12. 
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FIG. 9: (color online) Plot of r + a/I5 vs. \r\avg for c = 0.80. The values of the exponent are extracted from the best collapse of 
the avalanche size distribution for each set of three a values ranging from a = 1.00 — 1.44 and different values of L = 168, 240, 390 
and 600. The extrapolated value is the point at \r\avg —t 0. 



FIG. 10: (color online) Variation of the exponent a with \r\avg based on the same data and procedure as used to obtain figure 
(9). 
























